!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!>\brief
!!
!! Define a state variable class.
!!
!! \n
!!
!! This module defines a state variable abstract type \c c_stv which
!! can be used as a building block for the implementation of linear
!! and nonlinear solvers, time integrators and so on. The main purpose
!! of this type is decoupling the details of the state variable
!! representation, such as internal fields and arrays with an
!! arbitrary number of dimensions, from the implementation of some
!! general purpose algorithms. Moreover, this allows a single \c c_sty
!! variable to be used in various such algorithms, avoiding copies.
!!
!! The operators defined for \c c_stv essentially are those required
!! by a vector space.
!<----------------------------------------------------------------------
module mod_state_vars

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_state_vars_constructor, &
   mod_state_vars_destructor,  &
   mod_state_vars_initialized, &
   c_stv

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members

 !> State variable
 !!
 !! \note As a general rule, the passed argument is the one which is
 !! set by the procedure.
 !! \note We avoid functions in the type fields because typically this
 !! would require frequent allocations/deallocations of the function
 !! results, while the subroutine arguments can be allocated once and
 !! then passed around thanks to the \c inout intent.
 !! \note In general, when extending the type \c c_stv allocatable
 !! and pointer fields will be added. In a general purpose subroutine,
 !! to make sure that all these fields are allocated correctly, local
 !! variables of type \c c_stv should be always allocatable and
 !! allocated with
 !! \code
 !!  allocate( local_var , source=input_arg )
 !! \endcode
 !!
 !! Note that for efficiency reasons it might be useful to override
 !! the implementation of the various operators provided here.
 !!
 !! We also need a comment concerning the parallel execution. In
 !! general terms, there are two kinds of distributed data: those for
 !! which shared data are assumed to be copies of the same value and
 !! those for which shared data are supposed to be added together. All
 !! the operations defined in \c c_stv are correct as far as all the
 !! operants are of the same type, which then is also the type of the
 !! result. For this reason, there is nothing in \c c_stv which refers
 !! explicitly to this aspect. However, it's the user's responsibility
 !! to avoid making operations with variables of different kinds, in
 !! which case the result is meaningless.
 type, abstract :: c_stv
 contains
  !> \f$x+=y\f$
  procedure(i_incr), deferred, pass(x) :: incr
  !> \f$x*=r\f$
  procedure(i_tims), deferred, pass(x) :: tims
  !> \f$x+=ry\f$ (requires a temporary)
  procedure,                   pass(x) :: inlt
  !> \f$z=x\f$
! compiler bug
! Workaround for the workaround for some older ifort versions
! The righ form is the following
!  procedure(i_copy), deferred, pass(z) :: copy
! The copy procedure should not be used: the abstract interface is
! better
  procedure, pass(z) :: copy
  !> \f$z=x+y\f$
  procedure,                   pass(z) :: add
  !> \f$z=rx\f$
  procedure,                   pass(z) :: mlt
  !> \f$z=x+ry\f$
  procedure,                   pass(z) :: alt
  !> \f$x=x+\sum_i r_iy_i\f$
  procedure,                   pass(x) :: inlv
  !> \f$z=\sum_i r_iy_i\f$
  procedure,                   pass(z) :: lcb
  !> Copy constructor (compiler bug)
  !!
  !! The contruct <code>allocate(y,source=x)</code> is still not
  !! completely supported: see DPD200180295 and
  !! http://software.intel.com/en-us/forums/showthread.php?t=103821
  !!
  !! The subroutine source is a temporary workaround where the copy is
  !! constructed by the user.
!  procedure(i_source), deferred, pass(x) :: source
! Workaround of the woraround: see the comments to i_copy for details
  procedure, pass(x) :: source
  procedure(i_source_vect), deferred, pass(x) :: source_vect
 end type c_stv

 abstract interface
  pure subroutine i_incr(x,y)
   import :: c_stv
   implicit none
   class(c_stv), intent(in)    :: y
   class(c_stv), intent(inout) :: x
  end subroutine i_incr
 end interface

 abstract interface
  pure subroutine i_tims(x,r)
   import :: wp, c_stv
   implicit none
   real(wp),     intent(in)    :: r
   class(c_stv), intent(inout) :: x
  end subroutine i_tims
 end interface

 abstract interface
  pure subroutine i_copy(z,x)
   import :: c_stv
   implicit none
   class(c_stv), intent(in)    :: x
   class(c_stv), intent(inout) :: z
  end subroutine i_copy
 end interface

 abstract interface
  pure subroutine i_source(y,x)
   import :: c_stv
   implicit none
   class(c_stv), intent(in)               :: x
   class(c_stv), allocatable, intent(out) :: y
  end subroutine i_source
 end interface
 abstract interface
  pure subroutine i_source_vect(y,x,m)
   import :: c_stv
   implicit none
   integer, intent(in) :: m
   class(c_stv), intent(in)  :: x
   class(c_stv), allocatable, intent(out) :: y(:)
  end subroutine i_source_vect
 end interface
 
 ! private members

! Module variables

 ! public members
 logical, protected ::               &
   mod_state_vars_initialized = .false.
 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_state_vars'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_state_vars_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
       (mod_kinds_initialized.eqv..false.   ) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_state_vars_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_state_vars_initialized = .true.
 end subroutine mod_state_vars_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_state_vars_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_state_vars_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_state_vars_initialized = .false.
 end subroutine mod_state_vars_destructor

!-----------------------------------------------------------------------

 pure subroutine inlt(x,r,y)
  real(wp),     intent(in) :: r
  class(c_stv), intent(in) :: y
  class(c_stv), intent(inout) :: x

  class(c_stv), allocatable :: tmp

   ! Compiler bug (ifort Version 12.1.1.256)
   ! The following two allocations should be equivalent, however when
   ! using the first one the input variable y gets modified by the
   ! call to mlt, possibly due to some poiter copy. The bug is
   ! DPD200180295 and is discussed here:
   ! http://software.intel.com/en-us/forums/showthread.php?t=103821
      !allocate(tmp,source=x)
      !allocate(tmp,source=y)
      call y%source(tmp)
   call tmp%mlt( r , y )
   call   x%incr( tmp )
   deallocate(tmp)

 end subroutine inlt
 
!-----------------------------------------------------------------------

 pure subroutine add(z,x,y)
  class(c_stv), intent(in) :: x, y
  class(c_stv), intent(inout) :: z

   call z%copy( x )
   call z%incr( y )

 end subroutine add
 
!-----------------------------------------------------------------------

 pure subroutine mlt(z,r,x)
  real(wp), intent(in) :: r
  class(c_stv), intent(in) :: x
  class(c_stv), intent(inout) :: z

   call z%copy( x )
   call z%tims( r )

 end subroutine mlt
 
!-----------------------------------------------------------------------

 pure subroutine alt(z,x,r,y)
  real(wp), intent(in) :: r
  class(c_stv), intent(in) :: x, y
  class(c_stv), intent(inout) :: z

   call z%mlt( r , y )
   call z%incr( x )

 end subroutine alt
 
!-----------------------------------------------------------------------

 pure subroutine inlv(x,r,y)
  real(wp), intent(in) :: r(:)
  class(c_stv), intent(in) :: y(:)
  class(c_stv), intent(inout) :: x

  integer :: i

   do i=1,size(r)
     call x%inlt( r(i) , y(i) )
   enddo

 end subroutine inlv

!-----------------------------------------------------------------------

 pure subroutine lcb(z,r,y)
  real(wp), intent(in) :: r(:)
  class(c_stv), intent(in) :: y(:)
  class(c_stv), intent(inout) :: z

  integer :: i

   call z%mlt( r(1) , y(1) )
   do i=2,size(r)
     call z%inlt( r(i) , y(i) )
   enddo

 end subroutine lcb

!-----------------------------------------------------------------------

! compiler bug
! The following subroutine should be eliminated and the abstract
! interface used instead.
  pure subroutine copy(z,x)
   class(c_stv), intent(in)    :: x
   class(c_stv), intent(inout) :: z

   ! This subroutine must be ALWAYS overridden!!
  end subroutine copy
  pure subroutine source(y,x)
   class(c_stv), intent(in)               :: x
   class(c_stv), allocatable, intent(out) :: y

   ! This subroutine must be ALWAYS overridden!!
  end subroutine source

end module mod_state_vars

